#this program determines which voxels lie within the boundary of the contour of interest
#those voxels within the contour boundary are marked by setting their dose value to zero 
import os #imports the operating system module, needed for file and directory functions
import numpy as np #imports the NumPy scientific computing library and makes np shorthand for numpy
os.chdir("/Users/Documents/pydicom-0.9.7") #sets the directory path
import dicom #imports the pydicom package
from matplotlib import path #imports the path module needed for determining if a point lies within a polygon
rtdosefile=dicom.read_file("composite dose.dcm") #reads in the dicom rtdose file
rtstructfile=dicom.read_file("structure set.dcm") #reads in the dicom rtstruct file
x_origin=float(rtdosefile.ImagePositionPatient[0]) #gets the dose array origin x coordinate
y_origin=float(rtdosefile.ImagePositionPatient[1]) #gets the dose array origin y coordinate
pixel_spacing=float(rtdosefile.PixelSpacing[0]) #gets the pixel spacing dimension
num_columns=rtdosefile.pixel_array.shape[2] #gets the number of columns in the dose array
num_rows=rtdosefile.pixel_array.shape[1] #gets the number of rows in the dose array
num_doseslices=rtdosefile.pixel_array.shape[0] #gets the number of slices in the dose array
num_contourslices=len(rtstructfile.ROIContours[0].Contours) #gets the number of contour slices
scaling_factor=0 #sets the scaling_factor to zero, this is used to mark the voxels that lie within the contour
contour=[] #initializes a list of arrays containing the coordinates of the polygon vertices for each contour slice
#the for loop below reads in an array containing the coordinates of the polygon vertices, for each contour slice
for contour_slicenumber in range(0,num_contourslices):
        contour.append(np.array(rtstructfile.ROIContours[0].Contours[contour_slicenumber].ContourData))
#the for loop below iterates through the slices of the dose array to determine which voxels lie within the contour
#boundaries.  Care must be taken to assure that the dose_slicenumber corresponds to the correct contour_slicenumber.
#In this particular case, dose_slicenumber = contour_slicenumber.
for dose_slicenumber in range(0,num_doseslices):
        #the line below deletes the z-coordinates of the polygon vertices
        xy_coordinates = np.delete(contour[dose_slicenumber], np.arange(2, contour[dose_slicenumber].size, 3))
        #the line below reshapes xy_coordinates into an Nx2 array of vertices needed by the path function
        vertices=np.reshape(xy_coordinates, (-1, 2))
        #the for loop below iterates through each voxel (aka pixel for a single slice) of a given slice in the dose
        #array and tests whether the voxel is within the contour boundary. If it is within the contour boundary, the
        #dose in the voxel is set to zero, thereby marking the voxels that lie within the contour.
        for col in range(0,num_columns):                   
                for row in range(0,num_rows):             
                        x=x_origin+col*pixel_spacing
                        y=y_origin+row*pixel_spacing
                        polygon = path.Path(vertices)
                        if polygon.contains_point([x,y]):
                                rtdosefile.pixel_array[dose_slicenumber,row,col]=    \
                                scaling_factor*rtdosefile.pixel_array[dose_slicenumber,row, col]
#the line below overwrites the PixelData attribute of rtdosefile
rtdosefile.PixelData=rtdosefile.pixel_array.tostring()
#the line below saves this file as "modified composite dose.dcm" which has 0 values for voxels within the contour
rtdosefile.save_as("modified composite dose.dcm")
rtdosefile_list=rtdosefile.pixel_array.tolist() #converts the rtdosefile array into a list
textfile=open("modified composite dose list.txt","w") #creates and opens a text file in write mode
print >>textfile, rtdosefile_list #writes the rtdosefile list to the text file
textfile.close() #closes the text file